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I .   INTRODUCTION 

A .   BACKGROUND 

Mankind  has  used  the  ideas  of  science  and  mathematics  to 
help  understand,  describe  and  control  the  world.  The  central 
notions  of  science  have  developed  gradually.  Truly  revolu- 
tionary ideas  do  not  often  surface.  The  development  of 
calculus  by  Newton  and  Leibnitz  in  the  seventeenth  century  is 
quite  possibly  the  last  great  revolutionary  development  in  the 
sciences.  Chaos  and  chaotic  dynamics  could  very  well  be  the 
next  scientific  revolution.  The  theories  of  chaotic  dynamics 
are  being  used  to  describe  and  interpret  many  natural  phenome- 
non that  have  heretofore  been  inexplicable. 

The  modern  study  of  chaotic  dynamical  systems  was  begun  by 
Henri  Poincare  at  the  beginning  of  the  twentieth  century. 
During  his  study  of  celestial  bodies,  Poincare  observed  that 
their  trajectories  were  forever  oscillating,  yet  irregular  and 
aperiodic.  Poincare  described  this  motion  as  chaotic.  He 
further  discovered  that  the  motion  of  a  chaotic  system  has  a 
very  sensitive  dependence  on  the  initial  conditions  of  the 
system.  Because  of  the  nonlinearity  of  chaotic  systems, 
further  rigorous  study  was  impossible  until  the  advent  of  the 
computer.  Since  1960,  many  mathematicians  and  engineers  have 
developed  techniques  to  analyze  chaotic  systems.  Much  debate 
has  occurred  about  the  nature  of  chaos.   Currently,  the  only 


universally  accepted  property  of  chaotic  dynamical  systems  is 
that  they  have  sensitive  dependence  to  initial  conditions. 
Agreement  about  certain  aspects  of  chaos  does  not  exist  within 
either  the  mathematical  community  or  the  engineering  community 
and  certainly  no  agreement  exists  between  these  groups. 

In  its  David  II  report  [Ref  1.1],  the  Committee  on  the 
Mathematical  Sciences  of  the  National  Research  Council,  listed 
chaotic  dynamics  as  one  of  the  brightest  areas  for  future  ad- 
vancement of  applied  mathematics.  The  report  states  that 
chaotic  behavior  has  been  found  in  nonlinear  systems  as 
diverse  as  ecology  and  economics  to  engineering  and  meteorolo- 
gy. The  theories  of  chaos  have  their  basis  in  the  mathemati- 
cal fields  of  topology,  measure  theory  and  combinatorics. 
Mathematicians  are  just  now  beginning  to  develop  numerical 
methods  for  computing  some  of  the  properties  which  Poincare 
described  some  eighty  years  ago. 

The  engineering  community  looks  at  chaos  from  a  different 
viewpoint .  Engineers  have  developed  techniques  to  analyze 
systems  that  have  discrete  data  points  to  determine  their 
chaotic  tendencies.  Pseudo  phase  planes,  Lyapunov  exponents 
and  multivariate  probability  scaling  analysis  are  but  a  few  of 
the  techniques  that  lend  insight  into  the  behavior  of  chaotic 
dynamical  systems.  No  single  method  is  completely  conclusive. 
Hence,  when  studying  nonlinear  systems,  a  variety  of  methods 
is  used. 


The  motivation  for  this  thesis  was  twofold.  First,  an 
attempt  has  been  made  to  bridge  the  gap  between  how  mathemati- 
cians develop  chaotic  analysis  and  how  engineers  apply  the 
methods  of  chaotic  analysis.  This  was  accomplished  by 
examining  systems  that  mathematicians  have  proven  to  be 
chaotic  and  using  engineering  methods  to  analyze  these 
systems.  The  fortran  program  CHAOS  developed  by 
CDR  Martinus  Sarigul-Kli jn  [Ref.  1.2]  while  studying  heli- 
copter vibrations  was  used  to  analyze  the  Lorenz  equations  and 
Duf f ing' s  equation.  Secondly,  the  data  generated  from  the 
Lorenz  equations,  Duf f ing' s  equation  and  the  helicopter 
vibration  data  used  by  CDR  Sarigul-Kli jn  was  analyzed  using  a 
new  technique  proposed  by  Osborne  and  Provenzale  [Ref  1.3] . 
All  computations  were  performed  on  a  VAX  based  system  in  the 
Advanced  Computational  Laboratory  of  the  Department  of 
Aeronautics  and  Astronautics,  Naval  Postgraduate  School. 

B.   OVERVIEW  OF  THESIS 

The  second  chapter  of  the  thesis  covers  background  needed 
to  understand  chaos.  The  chapter  begins  with  a  short  develop- 
ment of  the  mathematical  basis  for  chaos.  Probability  theory, 
specifically  the  central  limit  theorem,  is  also  developed. 
Then,  the  engineering  techniques  used  to  analyze  chaos  are 
presented.  Pseudo  phase  planes,  fractal  dimension,  Poincare' 
sections  and  Lyapunov  exponents  are  presented  here.  Finally, 
the  chapter  ends  with  a  literature  review  that  discusses 
recent   developments   and   relevant   papers    concerning 
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applications  of  chaos. 

The  third  chapter  begins  with  a  discussion  of  a  technique 
for  simulating  a  continuous  system  using  a  digital  computer. 
Next,  the  work  of  CDR  Sarigul-Kli jn  is  presented.  Next,  the 
methods  used  in  the  present  analysis  are  presented.  The 
multivariate  scaling  analysis  is  presented  here.  Further,  the 
method  used  to  analyze  a  continuous  system  using  a  discrete 
approximation  is  discussed. 

Chapter  IV  presents  a  summary  of  the  results  of  the 
analysis.  First,  the  Lorenz  equations  and  Duf f ing' s  equation 
are  evaluated  using  current  techniques.  Next,  the  multivari- 
ate scaling  analysis  results  are  presented  for  the  helicopter 
vibrations  data. 

Chapter  V  is  a  discussion  of  conclusions  from  the  analy- 
sis.  Scope  for  further  research  is  also  presented. 


II.   UNDERSTANDING  CHAOS 

A.  INTRODUCTION 

This  chapter  gives  a  brief  overview  of  the  background 
needed  to  understand  the  concept  of  chaos.  First,  the 
mathematical  basis  for  choatic  dynamical  systems  is  presented. 
Next,  the  engineering  techniques  used  to  analyze  chaos  are 
presented.  Finally,  the  chapter  ends  with  a  review  of  recent 
literature  that  concerns  applications  of  chaos  that  are 
relevant  to  the  present  study. 

B.  MATHEMATICAL  BASIS  FOR  CHAOS 

The  study  of  chaotic  systems  begins  with  the  study  of  dis- 
crete dynamical  systems.  Since  data  can  only  be  collected  at 
discrete  points  in  time,  it  is  quite  natural  to  examine  a 
dynamical  system  at  those  discrete  points  so  that  a  model  of 
the  system  can  be  constructed.  The  computer  has  allowed  the 
analysis  to  be  done  with  shorter  and  shorter  time  intervals. 
A  simple  example  of  a  discrete  dynamical  system  is  an  iterated 
function  system  constructed  with  a  calculator.  A  value  in 
this  system  depends  on  the  previous  value  in  the  system.  The 
system  could  be  constructed  by  choosing  a  value  and  then 
repeatedly  pressing  the  same  function  key  on  a  calculator.  If 
the  key  selected  were  the  square  root  function  and  the 
starting  number  was  100,  then  we  could  construct  a  sequence  of 


numbers  based  on  the  function  (square  root)  and  the  initial 
condition.  Clearly  this  system  is  quite  predictable  but  it  is 
still  dependent  on  the  initial  conditions.  A  system  with  this 
property  is  said  to  be  deterministic. 

Dynamical  systems  can  be  classified  by  the  nature  of  the 
model  that  best  describes  that  system.  A  'linear'  system  is 
modeled  by  a  funtion  that  is  a  straight  line  that  passes 
through  origin.  A  function  is  called  'affine'  if  its  graph  is 
a  straight  line  that  does  not  pass  through  the  origin.  An 
affine  system  could  be  modeled  by  the  relationship: 

Xn+1  =  AXn  +  B. 
A  'nonlinear  system'  takes  the  form: 

Xn+1  =  f  (Xn)+g(n) 
where  f (Xn)  =  AXn  and  g(n)  is  a  function  of  the  iteration 
number.  The  graph  of  a. nonlinear  system  is  not  a  straight 
line  and  need  not  pass  through  the  origin.  For  example,  if 
f(X)  =  x  and  g(n)  =  n  the  equation  for  the  15th  iteration 
would  be 

X15=f  (X14)=X14  +15. 
A  simple  iterated  function  system  might  revisit  a  point 
after  a  certain  interval.  If  Xn+1  =  f(Xn)  =  Xn  then  Xn  is 
called  a  'fixed  point'.  If  Xn  =  f  (Xn)  =  Xn+m  then  Xn  is  called 
a  'periodic  point'.  That  is,  if  one  starts  at  point  Xn  and 
performs  'm'  iterations  on  the  system  the  last  iteration  will 
yield  point  Xn.  A  fixed  point  remains  fixed  after  successive 


iterations,  while  a  periodic  point  is  revisited  only  after  a 
certain  number  of  function  iterations. 

Periodic  points  and  fixed  points  can  be  classified  as 
attracting  or  repelling.  A  point  is  called  'repelling'  if  a 
small  perturbation  in  the  initial  condition  causes  the 
function  to  diverge  exponentially.  A  point  X  is  called  at- 
tracting if  a  small  perturbation  in  the  initial  condition 
causes  the  system  to  converge  to  point  X  after  successive 
iterations.  That  is,  if  a  system  is  iterated  with  two  sets  of 
initial  conditions,  one  condition  an  attracting  point  and  one 
condition  £  distant  from  the  attracting  point,  then  after 
successive  iterations  the  functional  values  will  be  closer 
together  than  the  original  points.  Figure  2.1a  and  2.1b 
exhibit  a  repelling  point  and  an  attracting  point,  respective- 
ly. It  can  be  shown  that  if  |  f '  (Xn)  |  <  1  then  Xn  is  an  attract- 
ing point.  While,  if  |  f '  (Xn)  |  >  1  then  Xn  is  a  repelling 
point . 

Before  giving  a  somewhat  generally  accepted  mathematical 
definition  of  chaos,  the  concepts  of  transitivity,  denseness 
and  other  basic  mathematical  concepts  must  be  explored.  Many 
good  references  are  available.  Devaney  [Ref.  2.1]  presents  a 
particularly  good  explanation  of  these  mathematical  concepts. 

First  the  definition  of  metric  space  is  presented.  A 
'metric  space'  (X,d)  is  a  space  X,  together  with  a  real  valued 
function  d,  that  maps  points  in  X  to  points  in  X,  and  which 


measures  the  distance  between  pairs  of  points  in  X.    The 
function  d  must  obey  the  following  rules: 

(1)  d(x,y)  =  d(y,x)  for  all  x  and  y  in  X 

(2)  0  <  d(x,y)  <  oo  for  all  x,  y  in  X,  x  *  y 

(3)  d(x,x)  =  0  for  all  x  in  X 

(4)  d(x,y)  <  d(x,z)  +  d(z,y)  for  all  x,  y,  and  z  in  X 

A  point  is  called  a  'limit  point'  if  there  is  a  sequence 
of  points  that  converge  to  that  point.  A  connection  between 
attractors  and  limit  points  becomes  evident.  The  concept  of 
convergence  is  common  to  both  attractors  and  limit  points. 
Next,  the  'closure'  of  a  subset  of  a  metric  space  is  said  to 
be  the  union  of  that  subset  and  all  of  its  limit  points. 
Next,  'denseness'  is  defined.  Let  (X,d)  be  a  metric  space. 
Let  B  be  a  subset  of  X.  B  is  said  to  be  dense  in  X  if  the 
closure  of  B  equals  X. 

A  set  is  'closed'  if  it  contains  all  of  its  limit  points. 
If  a  set  is  not  closed  it  is  said  to  be  'open' .  The  concept 
of  'transitivity'  is  next  presented.  Let  U  and  V  be  open 
subsets  of  the  metric  space  (X,d) .  A  dynamical  system  {x,f} 
in  that  metric  space  is  transitive  if  there  exists  a  finite 
integer  such  that 

U  n  fon(V)  *0 

In  simple  terms,   this  means  that  any  arbitrarily  small 
neighborhood  in  X  will,  under  successive  iterations,  overlap 
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with  any  other  arbitrarily  small  neighborhood  in  X.  There- 
fore, the  system  cannot  be  divided  into  two  disjoint  open 
sets . 

Now  the  definition  of  'chaos'  as  defined  in  [Ref.  2.1]  is 
presented.   A  system   is  said  to  be  chaotic  if: 

(1)  it  has  sensitive  dependence  on  initial  conditions 

(2)  it  is  topologically  transitive 

(3)  periodic  points  are  dense. 

An  important  phenomena  in  the  study  of  chaos  is  the 
'bifurcation'.  In  general,  dynamical  systems  arrive  at  a 
chaotic  state  through  a  bifurcation  of  one  type  or  another. 
Bifurcations  are  found  throughout  the  real  world.  A  bifurca- 
tion is  simply  a  dramatic  change  in  a  dynamical  system.  Water 
freezing,  beam  buckling  under  loads  and  the  splitting  of  an 
atom  during  a  nuclear  explosion  are  all  examples  of  bifurca- 
tions. Many  systems  are  nonchaotic  until  a  parameter  or  set 
of  parameters  attain  certain  critical  values.  As  the  system 
achieves  a  critical  state  a  bifurcation  occurs  and  the  system 
becomes  chaotic.  In  general,  not  all  bifurcations  are  identi- 
cal . 

A  'saddle-node  bifurcation'  occurs  when  a  function  has  a 
semi-stable  fixed  point.  A  semi-stable  fixed  point  is 
attracting  if  perturbed  in  one  direction  and  repelling  if 
perturbed  in  another  direction.  In  the  simplest  two  dimen- 
sional case,  the  bifurcation  occurs  as  the  function  approaches 


the  line  y  =x.  When  the  graph  of  the  function  just  touches 
the  line  the  function  splits  into  two.  Figures  2.2a  and  2.2b 
provide  graphical  representations  of  a  saddle-node  bifurcation 
of  the  logistic  equation. 

A  period  doubling  bifurcation  occurs  when  an  attracting 
fixed  point  becomes  repelling.  As  the  parameter  achieves  its 
critical  value,  the  fixed  point  now  becomes  a  cycle  of  period 
two.  As  the  parameter  passes  another  critical  value,  yet 
another  split  occurs.  This  process  continues  until  a  chaotic 
state  is  achieved.  Saddle-node  and  period  doubling  are  the 
two  most  common  bifurcations.  While  others  are  found  in 
nature,  it  is  clear  that  a  route  to  chaos  is  through 
bifurcations  [Ref.  2.2]. 

A  remarkable  property  about  dynamical  systems  was  proved 
by  Sarkovskii  [Ref  2.1].  Sarkovskii  developed  a  procedure  to 
determine  the  periods  of  a  dynamical  system.  He  showed  that 
by  ordering  the  integers  in  a  certain  manner,  the  periodic 
points  of  a  function  could  be  predicted.  As  previously 
discussed,  a  periodic  point  is  a  point  that  is  revisited  after 
a  number  of  functions  iterations.  In  this  ordering  if  a 
system  has  a  periodic  point  of  a  higher  order  period  then  it 
has  periodic  points  of  all  the  periods  that  follow.  For 
example,  if  a  system  has  a  point  of  period  three,  then  it  has 
periodic  points  of  all  periods.  To  find  points  of  other 
periods  one  must  change  the  initial  conditions.  The  highest 
order  period  is  three  and  the  lowest  order  period  is  one. 
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Sarkovskii' s  theorem  states   as  follows: 

let  the  ordering  of  the  integers  representing  periodic 
points  of  a  system  be: 

3>5>7>9 ...>... 
2x3>2x5>2x7> ...>... 
22x3>22x5>  ...>... 
25>2*>23>22>2>1 


then,  if  a  function,  f,  has  a  periodic  point  of  period  n,  and 
n  >  K  in  the  ordering,  then  f  also  has  a  periodic  point  of 
period  K.  Hence,  period  3  implies  all  other  periods.  Period 
2  only  implies  period  1 . 

The  next  concept  to  be  discussed  in  this  section,  is  the 
notion  of  a  fractal.  Fractals  are  subsets  of  metric  spaces. 
Fractals  have  become  important  in  helping  to  classify  and 
quantify  chaos.  Fractal  dimension  is  a  useful  tool  in 
measuring  chaos. 

The  mathematics  behind  the  fractal  is  quite  rigor- 
ous. [Ref .2.3]  A  basic  example  will  present  a  clear  picture  of 
a  fractal.  The  Cantor  middle  thirds  set  is  the  archetypical 
fractal.  This  set  is  constructed  by  starting  with  the 
interval  [0,1]  .  The  function  removes  the  middle  third  of  each 
remaining  segment.  After  successive  iterations  the  fractal 
remains.  Figure  2.3  illustrates  the  Cantor  set.  What  remains 
is  the  attractor  for  the  system.  Each  segment  is  similar  to 
the  previous  segment.  This  important  property  is  called  self- 
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similarity.  The  function  for  generating  a  fractal  can  be  a 
rotation,  a  translation,  a  reflection  or  any  combination  of 
those  three.  Many  functions  or  mappings  may  be  needed  to 
generate  a  particular  fractal.  In  all  cases,  the  function  (or 
functions)  acts  on  the  previous  point  or  segment.  However,  to 
generate  a  fractal  each  function  must  be  contractive  rather 
than  expansive.  This  means  that  the  function  cannot  enlarge 
the  previous  segment.  But,  the  combination  of  multiple 
functions  can  be  larger  than  the  original  segment. 

'Fractal  dimension'  provides  a  useful  measure  of  a  chaotic 
system.  A  chaotic  system  will  have  a  noninteger  fractal 
dimension.  While  not  all  systems  with  noninteger  fractal 
dimensions  are  chaotic,  this  is  still  a  useful  measure  [Ref. 
1.3]  . 

The  fractal  dimension  is  a  measure  of  how  much  space  an 

object  occupies  in  its  metric  space.   An  e-ball  about  a  point 

is  simply  all  points  less  than  £  distance  from  the  point.   To 

cover  a  set,  one  should  take  the  smallest  number  of  e-balls 

that  intersect  all  points  in  the  set.   The  fractal  dimension 

is  defined  as  follows: 

Let  A  be  an  element  of  the  Hausdorf  space  where  (X,d)  is 
the  metric.  For  each  £  >  0  let  N (A,  E)  denote  the  smallest 
number  of  closed  balls  of  radius  £  >  0  needed  to  cover  A. 

D  _  Lim       Ln(N(A,e)  ) 


e^O     Ln(l/e) 
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Then,  element  A  is  said  to  have  fractal  dimension  D.  The 
fractal  dimension  gives  a  measure  of  the  relative  size  of  the 
attractor  of  a  chaotic  system. 

The  examples  discussed  above  were  all  one  dimensional  in 
nature.  However,  the  concepts  of  chaotic  dynamics  can  easily 
be  applied  to  systems  of  higher  dimensions. 

The  last  concepts  discussed  in  this  section  are  the 
relevant  areas  of  probability  theory.  The  probability  density 
function  is  the  relative  measure  of  the  frequency  of  occur- 
rences. For  a  discrete  random  variable  the  summation  of  all 
values  of  the  probability  density  function  is  unity.  That  is, 
if  f (xn)  is  the  probability  of  occurrence  of  point  xn  then 

N 

72  =  1 

This  concept  is  intuitively  obvious.  The  analog  to  summation 
for  a  continuous  system  is  the  integral;  for  a  system  with 
more  than  one  random  variable  the  concept  does  not  change. 
For  a  multivariable  system  the  summation  over  all  random 
variables  will  equal  unity.  The  summation  over  a  single 
random  variable  of  a  multivariable  system  will  always  be  less 
than  or  equal  to  unity. 

Data  can  be  distributed  in  a  variety  of  patterns.  If  each 
value  in  the  range  has  an  equal  probability  of  occurring,  the 
data  is  said  to  have  a  uniform  distribution.  If  the  data  is 
symmetrically  dispersed  about  a  central  point  with  more  data 
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near  the  center  of  the  range  than  at  the  ends  of  the  range, 

the  data  has  a  Gaussian  or  normal  distribution.   Figure  2.5 

illustrates  some  common  probability  distributions.    A  key 

theorem  allows  data  with  any  distribution  to  be  transformed  to 

normal  distribution.   The  Central  Limit  Theorem  follows: 

Theorem:  Let  X:,  .  .  .  Xn  be  independent  and  identically 
distributed  random  variables  all  with  expectations  =  |!  and 
variance  equal  to  a 


y   =  ^^ 


then  Yn  converges  to  a  standard  normal  random  variable. 

This  theorem  is  applied  by  Osborne  and  Provenzale  [Ref  1.3] 

when   comparing   chaotic   systems   with   randomly  generated 
systems . 

C.   ENGINEERING  TECHNIQUES 

Most  naturally  occurring  phenomenon  have  a  continuous 
nonlinear  pattern.  To  study  these  events,  engineers  must  take 
discrete  readings  over  a  time  interval  or  across  a  distance. 
The  very  act  of  measurement  applies  a  strain  on  the  system  and 
slightly  changes  the  system.  The  process  of  discretizing  the 
continuous  data  produces  an  imperfect  model  of  the  system. 
Again  error  is  introduced  into  the  data.  This  section  will 
discuss  the  engineering  techniques  used  in  the  study  of 
chaotic  dynamics. 


14 


1 .  Time  Domain  Analysis 

The  three  steps  to  convert  an  analog  signal  to  a 
digital  signal  are  sampling,  quantizing  and  encoding  [Refs. 
1.2  and  2.4],  Sampling  a  continuous  signal  produces  a  series 
of  discrete  values.  A  generally  accepted  practice  is  to 
sample  the  data  at  twice  the  rate  of  the  uppermost  frequency 
of  the  signal.  The  sampled  data  must  be  quantized.  The 
spectrum  is  divided  into  different  levels  known  as  quantum 
levels.  The  signal  is  then  compared  with  each  quantum  level 
and  the  reading  at  a  specific  time  is  recorded  at  the  nearest 
quantum  level.  The  maximum  error  introduced  here  is  then  one- 
half  the  quantum  level.  The  next  step  is  translating  the 
quantum  levels  to  the  maximum  number  of  levels  available  in 
the  data  recording  equipment.  This  is  know  as  machine 
precision.   Again  an  error  is  introduced  into  the  data. 

2.  Frequency  Domain  Analysis 

An  alternative  method  for  analyzing  a  continuous  data 
signal  is  a  frequency  domain  analysis.  This  method  uses  a 
Fourier  analysis  to  measure  a  function  using  harmonic  compo- 
nents that  have  varying  amplitudes,  frequencies  and  phases. 
A  signal  in  the  frequency  domain  can  be  represented  by  two 
plots.  A  power  spectral  density  graph  is  a  plot  that  is  often 
used  to  represent  signals.  The  square  of  the  amplitude  is 
plotted  against  the  frequency. 
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3.   Phase  Planes  and  Pseudo  Phase  Planes 

The  concept  of  the  phase  plane  was  introduced  by  Poin- 
care.  A  phase  plane  represents  a  dynamical  system  by  using 
two  independent  properties  that  best  describe  the  system 
dynamics.  At  each  point  in  time  those  properties  are  measured 
and  then  plotted  on  a  two  dimensional  graph.  At  any  point  in 
time  the  system  can  be  described  by  a  single  point  on  the 
graph.  The  resulting  plot  is  used  to  analyze  the  system.  If 
the  system  approaches  a  limit  point  or  limit  cycle  over  time 
it  is  said  to  have  an  attractor.  The  phase  plane  requires  two 
separate  measurements  of  a  system.  As  was  discussed  earlier, 
each  measurement  introduces  an  error  into  the  data.  Also, 
collecting  data  for  two  variables  will  be  more  costly  and 
might  be  physically  impossible.  Hence,  a  phase  plane  analysis 
may  be  impractical. 

Takens  suggested  a  simpler  technique  for  constructing  a 
phase  plane  diagram  [Ref.  2.5].  This  technique  is  called  a 
pseudo  phase  plane.  A  dynamical  system  is  measured  with  only 
one  variable.  This  procedure  generates  a  single  time  series. 
A  second  time  series  is  constructed  by  adding  an  embedding 
time  to  the  original  data.  If  the  embedding  time  is  two  time 
periods,  the  third  point  in  the  original  series  becomes  the 
first  point  in  the  new  series.  The  data  in  the  second  series 
is  then  plotted  versus  the  first  series  to  construct  a  pseudo 
phase  plane.    The  pseudo  phase  plane  captures  the  system 
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dynamics  without  filtering  the  signal  by  integration  or 
differentiation  of  the  measured  signals  [Ref.  2.6].  Any 
dynamic  property  of  the  system  can  be  used  to  represent  the 
system.  The  choice  of  embedding  values  plays  a  critical  role 
in  the  accuracy  and  interpretation  of  the  analysis.  A  small 
embedding  time  produces  a  centralized  data  set  while  a  large 
embedding  time  introduces  error  through  edge  effects. 
Grassberger  [Ref. 2. 7]  suggests  that  the  embedding  time  should 
be  about  one  quarter  of  the  most  predominant  frequency  of  the 
data.  Dvorak  and  Klaschka  [Ref.  2.8]  have  since  proven  that 
an  embedding  time  between  five  and  eleven  time  intervals 
produces  a  pseudo  phase  plane  with  the  least  amount  of  bias. 
In  recent  years,  many  methods  have  been  developed  for 
reducing  the  measurement  error  or  noise  in  the  analysis. 
Hammel  [Ref.  2.9]  has  been  able  to  reduce  the  effects  of  noise 
tenfold.  This  produces  a  data  set  that  is  accurate  to  machine 
precision  on  even  the  most  advanced  computers.  For  an 
experimental  procedure  or  a  numerical  simulation,  the  pseudo 
phase  plane  is  more  efficient,  less  costly  and  produces 
results  that  are  just  as  accurate  as  the  phase  plane  and 
captures  the  original  system  dynamics. 

Pseudo  phase  planes  can  be  constructed  in  higher  dimen- 
sions as  well.  The  procedure  to  construct  higher  dimensional 
pseudo  phase  planes  is  similar  to  that  used  to  construct  the 
two  dimensional  phase  planes.   The  third  and  subsequent  time 
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series  are  constructed  using  the  original  time  series  and  new 
embedding  factors.  The  embedding  factor  should  be  different 
for  each  new  time  series.  More  than  three  dimensions  cannot 
be  depicted  graphically  but,  multidimensional  pseudo  phase 
planes  can  provide  insight  into  the  dynamics  of  a  system. 

Poincare  sections  can  be  constructed  from  these  phase 
planes  [Ref.  2.10].  A  Poincare  section  is  a  two  dimensional 
cut  of  a  three  dimensional  phase  plane.  The  section  shows 
where  trajectories  of  the  system  pierce  the  phase  plane.  This 
section  allows  one  to  graphically  determine  if  a  bounded 
attractor  or  limit  cycle  is  present.  A  well  defined  attractor 
is  evidence  of  a  chaotic  system.  The  Poincare  section  can  be 
taken  at  any  angle  so  that  the  number  of  possible  sections  is 
infinite.  Higher  dimensional  Poincare  sections  may  be 
obtained  similarly. 

4 .   Van  Der  Pol  Planes 

In  higher  dimensional  space,  Poincare  sections  can  be 
obtained  as  intersections  of  the  trajectories  with  a  specified 
hyper-plane.  The  internal  structure  of  the  attractor  may  then 
be  examined.  The  Van  Der  Pol  plot  is  constructed  by  untwisting 
the  trajectory  at  a  certain  rate  to  create  a  series  of 
Poincare  sections.  The  sections  are  then  merged  together  to 
yield  a  Van  Der  Pol  plot  [Ref.  2.11].  These  are  particularly 
useful  when  a  single  frequency  dominates  the  response  of  a 
system.   More  information  on  Poincare  sections  can  be  found 
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in  [Ref .  2.6]. 

5.   Lyapunov  Exponents 

The  methods  discussed  previously  have  all  provided  a 
graphical  representation  of  the  chaotic  nature  of  dynamical 
systems.  The  graphs  offer  pictorial  evidence  of  the  chaotic 
nature  of  the  system  being  analyzed.  These  methods  do  not 
provide  quantitative  measures  of  chaos.  The  Lyapunov  exponent 
offers  quantitative  evidence  of  chaos.  Wolf  [Ref.  2.12] 
offers  an  algorithm  for  calculating  these  exponents  based  on 
a  time  series.  This  method  measures  the  attraction  or  repul- 
sion over  time  of  two  adjacent  trajectories  with  different 
initial  conditions.  A  positive  exponent  means  that  the 
distance  between  the  systems  has  grown  greater  with  time.  A 
negative  exponent  indicates  that  the  trajectories  have  grown 
closer.  A  positive  exponent  is  an  indication  of  unpredict- 
ability and  chaotic  behavior  in  a  dynamical  system.  A 
positive  exponent  alone  is  not  enough  to  determine  if  a  system 
is  chaotic.  Recent  studies  by  Kapitaniak  and  Naschie  [Ref. 
2.13]  show  that  a  purely  random  system  can  also  have  positive 
Lyapunov  exponents. 

The  Lyapunov  exponents  also  gives  a  measure  of  the  rate  of 
information  loss  in  a  system.  As  such,  the  exponent  can  be  a 
useful  gage  as  to  how  far  into  the  future  one  can  predict  the 
dynamics  of  a  system. 
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6.  Fractal  Dimension 

Another  measure  used  to  quantify  chaos  in  a  dynamical 
system  is  the  fractal  correlation  dimension.  As  discussed 
earlier,  fractal  dimension  gives  a  lower  bound  for  the  number 
of  degrees  of  freedom  needed  to  model  a  physical  system.  A 
finite  non-integer  fractal  dimension  suggests  yet  further 
evidence  of  the  presence  of  a  strange  attractor  and  chaos. 
However,  Osborne  and  Provenzale  have  shown  that  a  randomly 
constructed  data  can  also  have  a  non-integer  fractal  dimension 
[Ref.  1.3].  Ruelle  [Ref.  2.14]  has  shown  that  all  fractal 
dimension  estimates  that  approach  21ogN  are  suspect.  In  this 
case,  N  is  the  number  of  data  points  the  time  series. 
Therefore,  for  most  time  series  of  10,000  data  points  fractal 
dimension  approaching  six  yields  inconclusive  evidence  of 
chaos  in  the  dynamical  system. 

As  with  much  of  this  field,  no  universally  accepted 
definition  of  fractal  dimension  exists  as  yet.  Farmer  et  al 
[Ref.  2.15]  describe  seven  natural  measures  of  dimension  that 
could  be  called  fractal  dimension.  These  measures  can  be 
divided  into  two  broad  classes.  One  class  requires  a  proba- 
bility measure  to  be  computed.  The  other  class  requires  only 
a  geometric  measurement. 

7 .  Evolving  Measures  of  Chaos 

Osborne  and  Provenzale  [Ref.  1.3]  have  suggested  that 
a  multivariate  scaling  analysis  will  lend  additional  insight 
into  the  dynamics  of  a  system.    A  multivariate  scaling 
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analysis  involves  computing  probability  density  functions  (PDF) 
of  the  time  series  data.  The  data  is  normalized  before  the 
PDF  is  constructed.  A  non-chaotic  system  will  generate  a 
Gaussian  PDF.  A  chaotic  system  will  yield  a  non-Gaussian  PDF. 
These  ideas  will  be  further  elaborated  in  Chapter  III. 

Lindsay  [Ref  2.16]  suggests  a  method  for  forecasting  data 
points  of  a  chaotic  system.  This  method  calls  for  calculating 
the  Jacobian  and  strange  attractor  to  use  global  rather  than 
local  properties  of  the  attractor.  A  predictor  function  is 
constructed  using  the  Jacobian  that  allows  one  to  forecast  as 
many  as  one-half  the  original  number  of  data  points. 

The  last  two  methods  utilize  the  global  properties  of 
chaotic  systems.  Locally,  random  systems  often  mirror  chaotic 
systems.  A  global  approach  helps  to  minimize  the  confusion 
caused  when  comparing  local  properties  of  random  systems  with 
local  properties  of  chaotic  systems. 

D.   SUMMARY 

Techniques  to  analyze  chaotic  dynamical  systems  have 
emerged  rapidly  during  the  last  decade.  Mathematicians, 
physicists  and  engineers  differ  in  their  approach  concerning 
analyzing  the  nature  of  chaotic  systems.  Methods  to  investi- 
gate chaos  in  systems  modeled  by  differential  equations  are 
well  documented.  Such  systems  are  continuous  and  need  not  be 
modeled  with  discrete  data  points.  Measuring  real  physical 
systems  yields  discrete  data  points  rather  than  continous 
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equations.  The  measurement  process  introduces  several  sources 
of  error  into  the  system.  Many  methods  can  help  describe  the 
chaotic  nature  of  such  a  system.  But,  as  a  result  of  the 
error  introduced,  no  method  gives  a  fool-  proof  description  of 
chaos.  Many  techniques  must  be  employed  concurrently  to 
ascertain  the  evidence  of  chaos  and  its  associated  properties 
that  allow  accurate  definition  of  the  dynamical  system. 

E.   LITERATURE  REVIEW 

Researchers  in  the  field  of  chaotic  dynamics  are  breaking 
new  ground  daily  and  technical  papers  concerning  chaos 
proliferate  in  leading  scientific  journals.  In  order  to  limit 
the  scope  of  this  literature  review,  the  papers  and  reports 
discussed  in  this  section  are  those  that  were  published 
between  December  of  1989  and  April  of  1991.  This  list  is  not 
exhaustive,  but  covers  the  papers  that  are  directly  related  to 
the  subject  matter  of  this  thesis. 

Bleher,  Grebogi  and  Ott  [Ref  2.17]  have  observed  a  new 
transition  from  a  system  in  a  steady-state  to  a  system  in  a 
chaotic  state.  This  transition  is  called  an  abrupt  bifurca- 
tion. The  abrupt  bifurcation  occurs  in  the  context  of  chaotic 
scattering.  The  behavior  is  caused  by  the  presence  of  a 
chaotic  invariant  set  of  unstable  bounded  orbits.  An  example 
showing  the  energy  involved  in  atom-diatom  collisions  is 
discussed.  Because  of  its  dependence  on  a  set  of  unstable 
bounded  orbits,  the  abrupt  bifurcation  might  help  explain  the 
nature  of  helicopter  vibrations. 
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Hughes  and  Proctor  [Ref .2.18]  studied  the  effects  of  noise 
on  a  chaotic  system.  They  show  that  noise,  even  at  very  low 
levels,  can  create  great  problems  in  determining  the  chaotic 
measures  of  a  system.  When  noise  was  added  to  a  known 
deterministic  chaotic  system  the  solution  to  the  system  is 
shown  to  be  indeterminate  and  incomputable. 

Corless,  Frank  and  Monroe  [Ref.  2.19]  determine  the 
chaotic  properties  of  the  Gauss  map.  The  Gauss  map  is 
developed  using  continued  fractions.  The  authors  determine 
the  fractal  dimension,  Lyapunov  exponent  and  the  correlation 
dimension  of  the  continous  Gauss  map.  All  of  these  properties 
indicate  that  the  Gauss  map  is  chaotic.  Next,  a  numerical 
computer  simulation  is  performed  to  generate  a  Gauss  map.  The 
simulated  map  has  fractal  dimension,  Lyapunov  exponent  and 
correlation  dimension  that  indicate  the  simulated  system  is 
nonchaotic.  The  simulated  system  is  created  to  machine 
precision.  The  observable  behavior  of  the  systems  is  not 
significantly  different.  This  suggests  that  a  system  with  a 
long  period  might  be  chaotic. 

A  numerical  experiment  using  a  self-affine  time  series  is 
presented  by  Higuchi  [Ref  2.20] .  The  randomness  of  the  phase 
distribution  strongly  affects  the  irregularity  of  the  system. 
This  conclusion  is  based  on  a  comparison  of  the  phase  distri- 
bution with  the  fractal  dimension.  The  randomness  of  the 
phase  distribution  is  shown  to  be  affected  by  the  behavior  of 
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the  system  in  the  time  domain. 

Liebert  and  Schuster  [Ref.  2.21]  report  the  effects  of 
varying  the  embedding  time  delay  when  analyzing  chaotic  time 
series.  For  an  infinite  time  series,  any  arbitrary  value  of 
an  embedding  value  is  acceptable  when  constructing  a  pseudo 
phase  plane.  However,  for  finite  time  series,  arbitrary 
choices  of  embedding  values  can  give  erroneous  results  when 
calculating  entropy  and  density.  A  method  for  calculating  a 
proper  embedding  value  which  requires  at  least  9000  data 
points  is  presented. 

Fang  [Ref.  2.22]  suggests  that  the  interplay  between  a 
large  noise  intensity  and  the  nonlinearity  of  the  system  may 
induce  a  nonequilibrium  transition  of  the  dynamical  behavior 
of  the  system.  The  degree  of  distortion  is  based  on  the  ratio 
of  the  noise  intensity  to  a  measure  of  the  system  nonlinear- 
ity. The  chaotic  measures  of  the  system  may  be  affected  by 
that  ratio.  Specifically,  the  Lyapunov  exponent  changes  sign 
as  the  noise  level  is  adjusted.  The  system  changes  from  a 
chaotic  state  to  a  periodic  state  with  the  introduction  of 
strong  white  noise. 

Kapitaniak  [Ref.  2.23]  presents  a  new  method  for  estimat- 
ing the  parameter  values  that  cause  a  dynamical  system  to  be 
chaotic.  This  method  applies  to  systems  that  exhibit  a  period 
doubling  route  to  chaos.  The  method  is  applied  to  Duffing' s 
oscillator.   The  method  is  based  on  an  approximate  analysis 
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and  the  Feigenbaum  universal  properties  of  period  doubling. 

Dvorak  and  Klaschka  [Ref  2.8]  present  a  modification  to 
the  Grassberger  and  Procaccia  algorithm  for  finding  correla- 
tion dimension.  The  original  algorithm  produces  a  bias  for 
embedding  dimensions  greater  than  five.  The  correlation 
dimension  of  the  system  is  changed  because  of  the  noise 
present  in  the  algorithm.  The  new  algorithm  allows  embedding 
dimensions  up  to  eleven  by  reducing  the  white  noise  level. 
This  algorithm  removes  the  bias  caused  by  the  edge  effects. 
All  other  causes  of  bias  from  the  old  algorithm  are  present  in 
the  new  algorithm. 

Chen  and  Ott  [Ref  2.24]  present  an  algorithm  for  computing 
the  cross-sections  of  strange  attractors.  These  cross- 
sections  are  useful  for  elucidating  the  geometry  of  the 
attractor.  The  cross-sections  also  allow  one  to  estimate  the 
fractal  dimension  of  higher  dimensional  attractors. 

Hammel  [2.9]  presents  a  method  for  reducing  the  noise  of 
a  time  series  of  discrete  data  points.  This  method  reduces 
the  noise  level  tenfold.  The  reduction  of  noise  allows  one  to 
distinguish  the  effects  of  noise  from  the  effects  of  chaos. 
Both  additive  noise  and  dynamic  noise  levels  are  reduced. 
This  method  produces  a  time-series  that  is  accurate  to  machine 
precision  and  thus  yields  more  accurate  calculations  of 
chaotic  measures  of  dynamical  systems. 

Stone  [Ref  2.25]  compares  the  effects  of  noise  versus 
periodic  forcing  on  the  power  spectra  of  a  Duffing  oscillator. 
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It  is  shown  that  the  form  of  the  power  spectral  density  of  the 
noise  driven  system  and  the  deterministic  chaotic  system  are 
indistinguishable.  This  suggests  that  a  high  frequency  power 
spectrum  cannot  distinguish  between  a  deterministic  and 
stochastic  system. 

Lindsay  [Ref.  2.16]  presents  a  simple,  computationally 
efficient  method  for  forecasting  chaotic  time  series  using  the 
Jacobian  and  certain  weighting  functions.  The  method  requires 
solving  a  small  number  of  linear  equations.  If  the  original 
time  series  has  10,000  data  points,  this  method  can  calculate 
the  next  5000  data  points  with  exceptional  accuracy.  The 
error  for  the  5000th  point  is  shown  to  be  less  than  2.7xl0"3%. 
The  same  method  can  be  used  for  noise  reduction. 
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Figure  2.1b   Attracting  Point 
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Logistic  Equation 


r    =    3.54 


Logistic  Equation 


r  =  3.56 


Figure  2.2  Period  doubling  route  to  chaos 
logistic  equation   F  (a) =r  (a)  (1-a) 
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Figure  2.3 

A  Typical  Fractal 
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III.   METHODOLOGY 

A.  INTRODUCTION 

This  chapter  presents  the  methods  used  to  generate  the 
simulated  data  and  the  methods  used  to  study  the  nonlinear 
dynamical  system  data  simulated  and  one  set  of  helicopter 
flight  vibration  data.  The  first  section  presents  the 
technique  used  to  generate  a  discrete  representation  of  a 
continuous  system  for  the  Lorenz  equations  and  Duffing' s 
equation.  A  discussion  of  randomly  generated  time  series  is 
also  presented.  Next,  Sarigul-Kli jn' s  work  is  discussed. 
This  section  focuses  on  the  calculation  of  the  fractal 
dimension  and  the  Lyapunov  exponent  of  a  dynamical  system 
represented  by  a  time  series.  A  development  of  the  multi- 
variate scaling  analysis  is  then  presented.  The  chapter 
closes  with  a  review  of  the  computer  simulations  used  in  the 
current  analysis. 

B.  SIMULATING  A  CONTINUOUS  SYSTEM 

In  1963,  Lorenz  with  the  help  of  a  computer,  studied  the 
atmosphere  of  a  convecting  fluid  modeled  by  the  now  famous 
Lorenz  equations . [Ref .  2.14]   These  equations  were  used 
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dx. 


dx7 

=px1-x2-x1x2 


dt 


where 

a  is  the  Prantl  number 

p  is  the  Rayleigh  number 

ft  is  the  aspect  ratio. 
Lorenz  studied  these  equations  with  the  following  combination 
of  parameters 

C  =  10 

p  =  28 

R  =  8/3 
used  to  model  the  system.  Lorenz  argued  that  the  atmosphere 
modeled  by  this  system  of  equations  is  chaotic  and  therefore, 
with  current  methods  longterm  behavior  could  not  be  predicted. 
The  earth's  atmosphere  is  similar  to  the  system  modeled  by  the 
Lorenz'  equations  and  thus,  longterm  weather  patterns  are  also 
unpredictable.  Lorenz  used  the  term  chaotic  to  mean  that  the 
system  had  sensitive  dependence  on  initial  conditions.  In 
1980,  Ueda  [Ref.  3.1]  studied  the  Duffing  equation  to  model 
the  motion  of  a  system  with  a  forcing  function  undergoing 
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large  elastic  deflection.   Ueda  found  that  the  time  series 
generated  from  the  behavior  of  this  equation 

*+0.05-^+x3=7.5cos(t) 


dt2  dt 

also  exhibited  chaotic  tendencies.  That  is,  a  long  term 
prediction  about  the  behavior  of  the  system  is  impossible. 
Both  Duf f ing' s  equation  and  the  Lorenz  equations  are  models 
of  continuous  dynamical  systems.  These  equations  also  model 
the  dynamics  of  several  systems. 

To  effectively  analyze  these  systems  with  a  digital 
computer  a  time  series  of  discrete  data  points  is  needed.  In 
the  present  study,  a  fourth  order  Runge-Kutta  numerical 
integration  was  used  to  generate  the  time  series.  Many 
excellent  references  are  available  for  this  and  other  numeri- 
cal integration  schemes  [Refs.  3.2  and  3.3]  .  The  fourth  order 
Runge-Kutta  technique  was  used  because  it  provided  the  level 
of  precision  required  for  the  analysis.  An  adaptive  step  size 
control  procedure  was  not  used  because  the  adaptive  strategy 
is  not  compatible  with  the  use  of  an  embedding  technique  for 
generating  a  pseudo  phase  plane.  To  accurately  generate  a 
pseudo  phase  plane  the  original  time  series  should  have 
uniformly  spaced  data.  Adaptive  step  size  does  not  provide 
evenly  spaced  data  and  therefore  was  impractical  for  the 
current  analysis.  Acceptable  precision  was  achieved  by  using 
a  step  size  of  .01  seconds.    The  equation  parameters,  step 
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size,  embedding  dimension  and  initial  conditions  were  all 
varied  to  generate  different  time  series  for  analysis.  The 
time  series  were  then  analyzed. 

Normally  distributed  colored  random  noise  was  numerically 
simulated  using  a  random  number  generator.  The  generator 
yielded  uniformly  distributed  random  numbers.  The  Central 
Limit  Theorem  was  applied  to  those  numbers  to  generate  a 
series  of  normally  distributed  numbers.  Twelve  uniform  random 
numbers  were  used  to  produce  one  Gaussian  random  number  [Ref . 
3.4]  .  The  numbers  were  then  added  to  a  sine  function  to 
produce  a  time  series  of  simulated  colored  noise.  The 
equation  used  to  generate  the  series  is 

X  =  .9+.25sin(2n) +.1RAND 

C.   AN  OVERVIEW  OF  SARIGUL-KLIJN' S  WORK 

Sarigul-Kli jn  performed  an  analysis  of  helicopter  flight 
vibration  data  [Ref.  1.2].  This  section  is  a  brief  descrip- 
tion of  part  of  his  work.  Although  helicopter  vibrations  are 
mostly  periodic,  evidence  of  chaos  was  claimed.  Specifically, 
the  Lyapunov  exponent  for  the  vertical  acceleration  measured 
under  the  pilot's  seat  was  found  to  be  0.3  to  1.7  bits  per 
second.  Also,  the  fractal  dimension  for  this  data  was  found 
to  converge  to  6.6  as  the  embedding  dimension  was  increased. 
The  fractal  dimension  for  the  lateral  and  longitudinal 
accelerations  converged  to  6.3  and  6.4,  respectively.  A 
random  signal  was  found  to  have  a  linearly  increasing  fractal 
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dimension.    This  random  signal  had  uniform  distribution, 
instead  of  the  more  widely  used  normal  distribution. 

The  Lyapunov  exponents  were  calculated  based  on  an  algori- 
thm developed  by  Wolf  et  al.  [Ref.  2.12].  The  exponent  is 
defined  as  follows 


«<■"»») 


where 

d(0)  is  the  initial  distance  between  trajectories 

d(t)  is  the  distance  at  time  't' 
The  algorithm  measures  the  distance  between  points  in  the 
original  time  series  and  the  points  in  the  embedded  time 
series.  A  positive  exponent  for  any  measurable  observable 
indicates  diverging  separation  between  trajectories  and  hence 
evidence  of  chaos . 

Sarigul-Kli jn  adopted  the  method  attributed  to  Grassberger 
and  Procaccia  to  calculate  fractal  correlation.  This  is  one 
of  many  measures  of  fractal  dimension  [Ref.  2.15].  This 
measure  of  fractal  dimension  is  defined  as  follows: 

C(r)=rd 
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where 


C(r)  is  the  probability  of  the  attractor  within  a  ball 

of  radius  r  (r-ball) 
r  is  the  radius  of  the  ball 
d  is  the  fractal  dimension. 
Therefore 


.  limf  InCir) 
a~r-0[      in  r 


The  fractal  correlation  dimension  is  computed  by  dividing  the 
number  of  points  inside  the  r-ball  by  the  total  number  of 
points  in  the  attractor.  The  process  is  repeated  for  several 
values  of  '  r'  .The  slope  of  the  In  C(r)  versus  In  r  curve  gives 
the  value  of  the  fractal  correlation  dimension. 

The  fractal  dimension  is  obtained  by  applying  this 
procedure  to  all  dimensions  of  the  pseudo  phase  space.  As  the 
phase  space  dimensions  increase,  the  correlation  dimension 
will  approach  an  asymptotic  value.  For  chaotic  systems,  this 
value  should  be  a  finite  noninteger  value. 


D.   MULTIVARIATE  SCALING  ANALYSIS 

This  section  presents  the  multivariate  scaling  analysis 
proposed  by  Osborne  and  Provenzale  [Ref.  1.3].  Osborne  and 
Provenzale  found  that  a  stochastic  process,  simple  colored 
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random  noise,  can  have  a  finite  fractal  correlation  dimension. 
Previously,  a  finite  fractal  correlation  dimension  was  thought 
to  be  positive  indication  of  the  presence  of  chaos.  A  new 
method,  the  multivariate  scaling  analysis,  was  developed  to 
aid  in  the  search  for  chaos. 

The  analysis  begins  with  a  time  series  of  data.  A  second 
time  series  is  generated  with  an  embedding  procedure  or  by 
taking  a  second  measurement  of  the  same  variable.  In  the 
present  analysis,  a  time  embedding  procedure  was  used  to 
generate  the  second  time  series.  The  two  time  series  are  used 
to  construct  a  scaling  space.  The  coordinates  for  the  scaling 
space  are  generated  as  follows 

^(t(i.1)M+J.)=xn(ti)-Xn(tJ),  for  i*j  (12) 

where 

Sn  are  the  scaling  space  coordinates 

Xn(ti)  is  the  original  time  series 

Xn(tj)  is  the  embedded  series 

M  is  the  number  of  points  in  the  original  series. 
Since  the  points  where  i^j  are  not  used,  the  time  series  of 
scales  will  have  M  -M  data  points.  The  embedding  dimension 
determines  how  many  sets  of  scaling  coordinates  are  generated. 
The  dimension  of  the  scaling  space  is  equal  to  the  embedding 
dimension.  If  the  embedding  dimension  is  two,  then  two  sets 
of  scaling  coordinates  are  generated.  A  second  series  of 
scaling  coordinates  is  constructed  using  the  same  procedure 
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but  a  different  embedding  value  than  the  first  series.  A 
joint  probability  density  function  is  then  constructed  using 
the  procedure  discussed  in  Chapter  II.  This  process  can  be 
applied  to  create  a  scaling  space  of  two  or  higher  dimensions. 
However,  a  graphic  analysis  only  allows  two  dimensions  to  be 
viewed.  The  probability  density  function  is  analyzed  to 
determine  the  chaotic  nature  of  the  system.  A  nonchaotic 
system  will  yield  a  joint  probability  density  function  that  is 
Gaussian  or  nearly  Gaussian  in  nature.  A  chaotic  system  will 
yield  a  joint  probability  density  function  that  is  clearly 
nonGaussian  in  nature. 

In  the  present  analysis,  multivariate  scaling  spaces  and 
joint  probability  density  functions  were  constucted  for  the 
Lorenz  equations,  Duffing' s  equation,  helicopter  vibration 
data  and  colored  random  noise.  While  the  multivariate  scaling 
analysis  does  not  provide  conclusive  proof  of  the  presence  of 
chaos,  it  does  provide  evidence  to  support  a  thorough  investi- 
gation of  a  dynamical  system. 

E.   COMPUTER  SIMULATIONS 

The  analysis  of  data  was  conducted  using  the  methods  de- 
scribed in  previous  sections.  Time  series  were  numerically 
generated  for  the  Lorenz  equations,  Duffing' s  equation  and 
colored  random  noise.  Many  different  parameter  values  were 
studied  for  the  Lorenz  and  Duf f ing' s  systems.  Different 
values  of  step  size  were  tried  for  the  numerical  integrations. 
For  all  cases,  effect  of  changes  in  embedding  value  was 
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studied.    Lastly,  for  all  cases,  various  sets  of  initial 
conditions  were  evaluated. 

The  evaluations  consisted  of  observing  time  history, 
computing  fractal  correlation  dimension,  Lyapunov  exponents, 
Fourier  spectrum  and  conducting  a  multivariate  scaling  space 
analysis  for  each  set  of  parameters,  embedding  values,  initial 
conditions  and  step  size.  Each  variable  could  conceivably 
affect  the  chaotic  nature  of  the  system.  Significant  results 
are  presented  in  Chapter  IV. 
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IV.   RESULTS 

A.   INTRODUCTION 

This  chapter  presents  the  results  of  the  analysis  of  the 
Lorenz  equations,  Duf f ing' s  equations,  random  colored  noise 
and  helicopter  flight  vibration  data.  The  results  are  first 
presented  for  each  dynamical  system.  In  order  to  study  the 
dynamical  system  described  by  a  system  of  ordinary 
differential  equations,  a  program  was  developed  based  on  a  4th 
order  Runge-Kutta  integration  scheme.  The  program  generated 
the  time  dependent  information  necessary  to  study  the 
dynamical  system.  Another  program,  based  on  the  multivariate 
scaling  analysis,  uses  the  time  histories  to  compute  joint 
probability  density  functions  as  discussed  in  earlier 
chapters.  Listings  of  these  programs  are  included  in 
Appendices  A  and  B,  respectively. 


B.   RANDOM  COLORED  NOISE 

Random  colored  noise  is  generated  using  the  following 
equation : 

X  =  .9  +  .25sin(27i)  +.1RAND 

The   data   generated   is   normalized   to   yield   a   Gaussian 
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distribution.  A  comparison  of  fractal  correlation  dimension 
for  random  colored  noise  and  purely  random  data  is  included  as 
Figure  4.1.  The  number  of  points  generated  is  10,000  with  an 
embedding  time  of  10  samples.  The  fractal  correlation 
dimension  of  the  random  data  increases  linearly  to  infinity 
with  the  embedding  dimension.  The  fractal  correlation 
dimension  of  the  random  colored  noise  increases  at  a  slower 
rate  than  the  embedding  dimension  and  appears  to  reach  an 
asymptotic  limit.  The  fractal  dimension  of  random  colored 
noise  is  a  noninteger  finite  number  for  all  values  of 
embedding  dimension  that  were  tested.  This  data  supports  the 
conclusion  of  Osborne  and  Provenzale  [Ref .  1.3] .  A  multivari- 
ate scaling  analysis  of  random  colored  noise  yields  a  joint 
probability  density  function  that  is  seen  to  be  Gaussian  in 
nature.   A  graph  of  this  function  is  included  as  Figure  4.2. 

C.   LORENZ  EQUATIONS 

The  time  series  for  the  Lorenz  equations  is  numerically 
generated  using  a  fourth  order  Runge-Kutta  integration 
numerical  scheme.  Step  sizes  of  .001  seconds,  .01  seconds  and 
0.1  seconds  were  used.  This  system  of  equations  could  not  be 
solved,  even  using  double  precision  arithmetic,  with  a  step- 
size  of  greater  than  .115  seconds  on  either  the  Micro-VAX  II 
computer  at  the  Naval  Postgraduate  School  or  the  CRAY-YMP 
computer  at  NASA-Ames  Center.   Step  sizes  larger  than  .115 
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caused  an  arithmetic  overflow  error.  Hence,  all  computations 
were  performed  with  step-sizes  less  than  .115  seconds.  Figure 
4.3  shows  a  comparison  of  the  X-Z  phase  plane  with  a  pseudo 
phase  plane  generated  using  the  X  time  series  only  and  an 
embedding  value  of  one.  10,000  data  points  were  used  to 
generate  the  pseudo  phase  plane.  The  pseudo  phase  plane  has 
the  same  general  structure  as  the  phase  plane  with  a  rotation 
of  one  branch  of  the  diagram.  Larger  embedding  dimensions 
caused  the  psuedo  phase  plane  to  rotate  and  generate  more 
branches.  The  fractal  correlation  dimension  increases  with 
step  size.  Figure  4.4  shows  a  comparison  of  fractal  dimension 
to  step  size.  A  step  size  of  .1  seconds  yields  a  fractal 
correlation  dimension  of  2.04  which  agrees  with  the  work  of 
Grassberger  and  Procaccia  [Ref.  4.1].  Figure  4.5  shows  plots 
of  the  slope  of  the  log(r)  versus  log  C(r)  curve  with 
differing  step  sizes. 

An  interesting  result  was  observed  as  the  Lorenz  equations 
were  studied.  Smaller  step  sizes  resulted  in  a  loss  of  low 
frequency  cycles  in  the  Fourier  spectrum.  Figure  4.6  shows 
the  Fourier  spectrum  plotted  to  the  Nyquist  frequency  for  the 
differing  step  sizes.  Figure  4.7  shows  plots  to  20%  of  the 
Nyquist  frequencies.  Figure  4.8  shows  the  low  frequency 
Fourier  spectrum  for  different  step  sizes.  The  smaller  step 
size  appears  to  act  as  a  filter  for  the  low  frequency 
components  of  the  Fourier  spectrum.  The  multivariate  scaling 
analysis,  using  a  10,000  data  points  and  a  time  delay  of  one 
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sample,  exhibited  a  non  Gaussian  behavior  suggesting  a  chaotic 
attractor  [Ref  1.3]  .  The  probability  distribution  function  is 
included  as  Figure  4.9. 

D.   DUFFING' S  EQUATIONS 

Next  Duffing' s  equation,  which  is  known  to  exhibit  chaotic 
behavior  is  studied.  Numerical  solution  was  obtained  by 
adapting  the  program  developed  for  the  Lorenz  equations. 
Figure  4.10  shows  a  comparison  of  the  X-Y  phase  plane  versus 
the  pseudo  phase  plane  generated  with  an  embedding  dimension 
of  1  applied  to  the  X  time  series.  Again,  fractal  correlation 
dimension  changed  with  step  size.  However,  in  this  case  the 
fractal  dimension  decreased  with  larger  step  sizes.  In  all 
cases,  the  fractal  correlation  dimensions  were  finite  non- 
integer  values.  The  range  of  fractal  dimensions  agrees  with 
the  work  of  Moon  and  Li'  [Ref.  4.2]  .  Figure  4.11  shows  a 
comparison  of  fractal  dimension  versus  change  in  step  size 
using  10,000  data  points  with  an  embedding  dimension  of  10  and 
a  time  delay  of  1  sample.  Figure  4.12  shows  the  Fourier 
spectrum  plotted  to  the  Nyquist  frequency  for  the  various  step 
sizes.  Figure  4.13  shows  the  Fourier  spectrum  plotted  to  20% 
of  the  Nyquist  frequency.  Figure  4.14  shows  a  comparison  of 
step  size  versus  the  low  frequency  Fourier  spectrum.  Again, 
the  smaller  step  size  acts  as  a  low  pass  filter.  But  in  this 
analysis  the  larger  step  size  appears  to  cause  the  data  to 
have  quasiperiodic  behavior.  Figure  4.15  shows  a  graph  of  the 
joint   probability   density   function   generated   from   the 
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mutlivariate  scaling  analysis.   The  graph  is  seen  to  be  non- 
Gaussian  indicating  chaotic  behavior. 

E.  HELICOPTER  VIBRATION  DATA 

This  section  presents  an  application  of  the  multivariate 
scaling  analysis  based  on  the  pseudo  phase  plane  method  for 
the  helicopter  vibration  data  [Ref .  1.2] .  The  vibration  data 
selected  for  the  present  study  are  the  vertical  acceleration 
recorded  under  the  right  pilot  seat  of  an  0H-6A  helicopter. 
Both  Higher  Harmonic  Control  'off  and  'on'  are  studied.  The 
joint  probability  density  function  with  a  time  delay  of  10 
data  points  is  included  as  figure  4.16.  The  graph  differs 
little  from  a  standard  Gaussian  distribution.  The  variation 
is  more  like  the  random  colored  noise  than  Lorenz  or  Duffing' s 
equations.  Figure  4.17  shows  the  joint  probability  density 
functions  for  all  of  the  dynamical  systems  studied.  This 
suggests  that  part  of  the  helicopter  vibrations  are  periodic 
while  part  of  the  vibration  is  chaotic.  This  conclusion 
supports  the  work  of  Sarigul-Kli jn  [Ref  1.2]. 

F.  CONCLUSIONS 

Two  trends  are  evident  in  the  analysis  of  the  data.  It 
appears  that  integration  step  size  filters  the  data  generated 
by  numerical  integration  using  a  computer.  Small  step  sizes 
might  cause  a  loss  of  the  low  frequency  components  while  large 
step  sizes  could  cause  the  behavior  of  a  chaotic  dynamical 
system  to  appear  to  be  quasiperiodic .    The  multivariate 
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scaling  analysis  appears  to  give  an  indication  of  the  chaotic 
nature  of  a  dynamical  system  distinguished  from  random  data. 
Still/  in  the  search  for  chaos,  this  technique  probably  does 
not  provide  a  definitive  answer  as  to  the  nature  of  chaotic 
behavior . 
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Figure  4.1   Fractal  dimension  vs.  embedding 
dimension  for  random  data  and  random  colored 
noise . 
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Figure  4.2   Joint  probability  density 
function  for  2-D  scaling  space  generated 
from  random  colored  noise  using  a  multi- 
variate scaling  analysis. 
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Figure  4.3a   X-Z  phase  plane  generated 
from  Lorenz  equations. 
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Figure  4 . 3b  X-x  pneudo  phase  plane 
generated  from  equations. 
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Figure  4.4  Fractal  dimension  vs. 
step  size 
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Figure  4.5  Slope  of  log(r)  vs. 
log  C(r)  for  Lorenz  equations. 
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dimension . 
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Figure  4.8  Low  frequency  Fourier 
spectrum  generated  from  Lorenz  equations 
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Figure  4.9   Scaling  space  joint 
probability  density  function  for 
the  Lorenz  equations. 
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Figure  4.10a   X-Y  phase  plane  for 
Duffing's  equation. 


Figure  4.10b   X-X  pseudo  phase  plane 
for  Duffing's  equation. 
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Figure  4.11  Slope  of  log(r)  vs.  log  C(r) 
for  Duffing's  equation.   The  maximum 
point  on  the  curve  is  the  fractal 
correlation  dimension. 
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Figure  4.12  Fourier  spectrum  plotted 
to  the  Nyquist  frequency  for 
Duffing' s  equation. 
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Figure  4.13  Fourier  spectrum  for 
Duffing' s  equation  plotted  to 
20%  of  the  Nyquist  frequency 
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Figure  4.14   Low  frequency  Fourier 
spectrum  for  Duffing's  equation. 


58 


Figure  4.15  Joint  probability 
density  function  for  scaling 
space  for  Duffing' s  equations 
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Figure  4.16  Joint  probability 
density  function  for  scaling 
space  for  helicopter  vibration 
vertical  acceleration  data. 
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Figure  4.17  Scalinq  space  joint 
probability  density  functions 
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V.   CONCLUSIONS  AND  SCOPE  FOR  FUTURE  RESEARCH 

Engineering  methods  of  chaos  theory  were  applied  to  a 
discrete  representation  of  chaotic  systems  of  equations.  The 
Lorenz  equations  and  Duf f ing' s  equation  were  examined  using  a 
time  series  approach.  The  time  series  were  generated  numeri- 
cally using  a  fourth  order  Runge-Kutta  integration  scheme. 
Several  simulations  were  executed  varying  the  equation 
parameters,  the  initial  conditions,  the  embedding  dimension 
and  the  integration  step  size.  Adaptive  step  size  control  was 
not  used  because  it  is  not  compatible  with  the  use  of  a  time 
embedding  procedure  to  generate  a  pseudo  phase  plane. 
However,  time  stepsize  affected  the  Fourier  spectrums  of  the 
systems  analyzed.  Small  step  size  acted  as  a  low  pass  filter 
in  the  chaotic  dynamical  systems. For  certain  parameters,  the 
systems  produced  Lyapunov  exponents  and  fractal  dimension  that 
indicated  the  presence  of  chaos.  The  size  of  the  embedding 
dimension  was  found  to  affect  the  chaotic  nature  of  the 
systems.  The  best  embedding  dimensions  were  found  to  be 
between  one  and  ten.  Large  embedding  dimensions  caused  the 
pseudo  phase  plane  to  have  a  different  nature  than  the  phase 
plane . 

A  new  method,  the  multivariate  scaling  analysis,  was 
applied  using  embedding  coordinates  to  the  Lorenz  equations 
and  Duf f ing' s  equation  as  well  as  to  the  helicopter  vibration 
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data  studied  by  Sarigul-Kli jn .  A  uniformly  distributed  random 
time  series  and  a  normally  distributed  random  time  series  were 
also  examined  using  this  new  technique.  When  examined  using 
the  multivariate  scaling  analysis,  chaotic  systems  yielded 
probability  density  functions  that  were  nearly  Gaussian.  The 
normally  distributed  random  time  series  produced  a  density 
function  that  was  Gaussian  in  nature.  Lastly,  the  helicopter 
vibration  data  yielded  a  density  function  that  was  non- 
Gaussian,  indicating  a  small  amount  of  chaos  in  helicopter 
vibrations . 

Many  avenues  are  open  to  further  research  in  this  area. 
Techniques  have  been  developed  that  allow  accurate  prediction 
of  future  points  in  a  chaotic  time  series.  This  technique  can 
be  used  to  predict  the  next  5000  points  in  a  time  series  for 
helicopter  vibrations.  This  technique  could  be  applied  to  the 
currently  available  data.  A  higher  harmonic  control  device 
utilizing  this  technique  could  be  used  to  actively  dampen  the 
helicopter  vibrations.  The  multivariate  scaling  analysis 
could  be  used  to  examine  other  experimental  data  that  appears 
to  have  a  chaotic  nature.  This  technique  offers  further 
evidence  of  the  chaotic  nature  of  the  system. 

The  Naval  Air  Test  Center  will  soon  instrument  an  H-60 
helicopter  with  a  higher  harmonic  control  system.  The  methods 
described  herein  could  be  applied  to  the  data  obtained  from 
that  study.   This  would  offer  further  evidence  of  the  chaotic 
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nature  of  helicopter  vibrations  and  help  to  determine  if 
vibration  data  can  be  accurately  predicted  and  hence  signifi- 
cantly reduced. 
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APPENDIX  A 
PROGRAM  TO  PERFORM  NUMERICAL  INTEGRATION 


*  This  program  solves  systems  of  ordinary  differential  equations 

*  using  a  fourth  order  Runge-Kutta  Scheme. 
* 

*  The  program  was  used  to  solve  both  Duffing's  equation  and  the 

*  Lorenz  equations. 
* 

* 

*  Two  subroutines  are  employed. 

*  DERIVS  calculates  the  values  for  the  ODE's 

*  RK4  is  the  Runge-Kutta  routine 
* 

*  As  written  the  program  does  not  use  an  adaptive  step-size 

* 

*  The  variables  are  defined  as  follows: 

*  N  -  the  number  of  ODE's  in  the  system 

*  NMAX  -  the  number  of  points  saved  for  plotting 

*  Currently  the  program  saves  every  tenth  point  in  a 

*  file  called  RKDAT.DAT 

*  H  -  the  desired  stepsize 

*  X,Y,Z  -  the  independent  variables 

*  Currently  the  program  solves  for  only  three 

*  variables.   If  a  system  has  more  than  three 

*  equations  the  additional  variables  must  be 

*  declared. 

*  F  -  the  ODE's 

*  the  equations  must  be  specified  in  the  subroutine 

*  '  DERIVS' 

*  Tim  -  time  stored  as  a  vector 
* 

* 
* 
* 

*  The  current  program  solves  the  duffing's  equation 

*  with  parameters: 

*  SIG,  R,  B 


All  other  variables  are  working  storage  spaces 


* 

*  Currently,  the  output  matches  the  HHC  data  from  Hughes 

* 
* 

*  The  Main  Program 

PARAMETER  (  N-=  2  ,  NMAX  =  9990  ) 

IMPLICIT  DOUBLE  PRECISION( A-H , 0-Z ) 

CHARACTER* 7  DUMMY 

DATA  H, A, B/. 1,7.5,0.5/ 

DIMENSION  XO(N) ,XEND(N) ,F(N),XKl(N),XK2(N),XK3{N),XK4(N),TO(10) 

COMMON  /PARS/H,A,B 

TIME=0.0 

********************************************  ********************* 

*  set  the  initial  conditions 
***************************************************************** 


DATA  XO/0. 1,4.1/ 
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DUM=U . U 
DUMMY='  AaaaaA' 

0PEN(8,FILE='dl' ,STATUS='old' ) 

**************************************************************** 
* 

*  the  lines  of  text  are  added  to  make  the  data  compatible 

*  with  the  program  'CHAOS' 
* 
*************************************************************** 

WR1TE( 8 , ' ( A) ' )  '  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx ' 
WRITE( 8 , ' ( A) ' )  '  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx ' 
WRITE( 8 , ' ( A) ' )  '  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx ' 

DO  11  I-1,NMAX 
X  =  XO( 1 ) 
Y=X0(2) 
TIME=TIME+H 

************This  statement  makes  the  data  compatible  with  ' CHAOS '*** ******** 

WRITE(8, 100)  DUMMY, TIME, X , Y , X , DUM , DUM 

CALL  RK4(XO,XEND, F , N , TIME , XKl , XK2 , XK3 , XK4 ) 

DO  14  L-1,N 

XO( L)=XEND(L) 
14  CONTINUE 

12  CONTINUE 

11     CONTINUE 
100     FORMAT (A7,F8 .3,Ell.4,4El2. 4) 

STOP 
END 

SUBROUTINE  RK4(XO,XEND, F , N , TIME , XKl , XK2 , XK3 , XK4 ) 

*  This  subroutine  performs  numerical  integration  using  4th  order 

*  Runge-Kutta  methods 

IMPLICIT  DOUBLE  PRECIS ION( A-H , O-Z ) 

COMMON  /PARS/H,A,B 

DIMENSION  F(N) ,XO(N) ,XEND(N) , XKl ( N ) , XK2 ( N ) , XK3 ( N ) , XK4 ( N ) 

*  DERIVS  returns  values  for  the  ODE'S 

*  Each  loop  is  a  step  for  N  equations 

CALL  DERIVS(XO,TIME,F,N) 
DO  51  K«1,N 

XKl(K)=H*F(K) 

XO(K)=XO(K)+XKl(K)/2.0D0 
51      CONTINUE 

CALL  DERIVS(XO,TIME,F,N) 
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XK2( K)=H*F( K) 

X0( K)=XO( K)+XK2( K)/2 . ODO 

52  CONTINUE 

CALL  DERIVS(XO,TIME, F,N) 
DO  53  K=1,N 

XK3(K)=H*F(K) 

XO( K)=XO( K)+XK3( K) 

53  CONTINUE 

CALL  DERIVS(XO,TIME,F,N) 
DO  54  K=1,N 

XK4(K)=H*F(K) 
XEND(K)=XEND( K)+(XKl(K)+2.0D0*XK2(K)+2 . 0D0 *XK3 ( K ) +XK4 ( K ) )/6 . 0D0 

54  continue 

RETURN 
END 


SUBROUTINE  DERI VS ( XO , T , F , N ) 

This  subroutine  obtains  values  for  the  ODE'S 

This  subroutine  must  be  changed  for  a  new  system  of  equations 

This  routine  solves  Duffing's  equation 

IMPLICIT  DOUBLE  PREC I SION ( A-H , O-Z ) 

COMMON  /PARS/H,A,B 

DIMENSION  F(N) ,XO(N) 

F(l)-XO(2) 

F(2)=A*COS(T)-XO( 1 )**3-B*XO(2) 

RETURN 
END 


*  use  this  routine  to  solve  the  Lorenz  equations 

******************************************************************** 

*  IMPLICIT  DOUBLE  PRECISION ( A-H , O-Z ) 

*  COMMON  /PARS/SIG,B,R,H 

*  DIMENSION  F(N),XO(N) 

*  F(l)«SIG*(XO(2)-XO(l) ) 

*  F(2)=R*XO(l)-XO(2)-XO(l)*XO( 3) 

*  F( 3)=XO(l)*XO(2)-B*XO(3) 

A 

*  RETURN 

*  END 

AAA******************************************************************* 
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APPENDIX  B 
MULTIVARIATE  SCALING  ANALYSIS 


This  program  performs  a  multivariate  scaling  analysis 

on  a  time  series  constructed  from  a  numerical  integration 

or  generated  from  experimental  results. The  program  uses  the 

central  limit  theorem  to  approximate  a  standard  gaussian 

distribution  so  that  if  the  data  is  truly  random  the  distribution 

will  be  random.   If  the  data  is  chaotic  the  data  will  be  skewed 

from  the  normal  distribution. 

parameters : 

ntot     total  values  in  sequence 

n       values  to  be  plotted.   fewer  values  are  plotted 

because  an  imbedding  process  is  used  to  generate  a 
sequence . 

j        the  number  of  divisions  on  each  axis 

x        the  initial  data  values 

xs       the  standardized  values 

xx       the  values  on  the  x-axis  for  the  graph 

yy       the  values  on  the  y-axis  for  the  graph 

s        sets  up  the  view  for  the  graph 

this  can  be  rotated(see  grafkit  manual) 

amp      the  function  values  for  the  graph  (z-values) 

work     a  workspace  grafkit  uses  to  generate  the  graph 
CHARACTERS  DUMMY 
CHARACTER*10  FNAME 

real  AMP (101, 101 ) , COUNT(101 ), Y( 3000), YS( 10000), TIM( 3000) 
real  X( 3000), XS( 10000), XX ( 101 ),YY( 101 ),S(6) ,WORK( 24440) 
DIMENSION  DUMMY ( 100),A(100),B(100),C(100),D(100) 
data  S/-8.0,-6. 0,3.0,0.0,0.0,0.0/ 

WRITE  (6,*)  'ENTER  THE  NUMBER  OF  DATA  POINTS' 
READ   (5,*)   NTOT 

WRITE  (6,*)  'ENTER  THE  EMBEDDING  VALUE' 

WRITE(6,*)   '    This  should  be  an  integer  value  between' 

WRITE(6,*)   '    0  and  100' 

READ  (5,*)  NSHIFT 

N=NTOT+NSHIFT 

variables : 

axis  the  discrete  values  on  the  axes 

1  the  number  used  for  the  embedding  process 

xs  the  original  data  values 

x  the  standardized  data  values 

y  the  embedded  data  values ( standard! zed ) 

istep  the  number  of  values  on  each  axis 


********************************************************************* 
*        read  a  file  here 

* 
********************************************************************* 
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WRITE(6,*)  'ENTER  THE  FILENAME' 

READ( 5, ' (A) ' )  FNAME 

0PEN(UNIT=8, FILE=FNAME,STATUS='OLD' ) 

DO  10  1=1, N 

READ(8,300)  tim(i),x(i) 

10  CONTINUE 

AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

A 
* 

*  As  written  this  format  statement  must  be  changed  to 

*  correspond  to  the  data  file  to  be  read 

A 
AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

300     FORMAT( t8, f 6 . 3, t25,el2.  4  ) 
CALL  OPNGKS( 1,2,1) 

AAAAAAAA***AAAAAAAAAAAAAAAAAAAAAAAAAAAA  +  AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

*  embedding  procedure 

AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA*AAAAAAAAAAA*A****AA*A*A*AAAAAAAA*AAAAAA* 

DO  11  I=l,NTOT 

Y( I )=X( I+NSHIFT) 

11  CONTINUE 

AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

A 

*  standardize  the  variables 

A 
AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

TOTAL=0.0 

DO  21  I=l,NTOT 

TOTAL~TOTAL+X( I ) 

21  CONTINUE 

RNTOT-NTOT 
XBAR=TOTAL/RNTOT 

SUM  =  0.0 

DO  2  2  I-l,NTOT 

SUM=SUM+( (X( I )-XBAR) **2) 

22  CONTINUE 

XSIGMA-SQRT( SUM/RNTOT) 

DO  23  I-l.NTOT 

X( I )=(X( I )-XBAR)/XSIGMA 

23  CONTINUE 

TOTAL=0.0 

DO  31  I«l,NTOT 

TOTAL=TOTAL+Y( I ) 
31       CONTINUE 
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SUM=0 . 0 

DO  32  1=1, NTOT 

SUM=SUM+( ( Y( I )-YBAR) **2 ) 

32  CONTINUE 

YSIGMA=SQRT( SUM/RNTOT) 

DO  33  I=l,NTOT 

Y( I )-( Y( I )-YBAR)/YSIGMA 

33  CONTINUE 

AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

*        construct  a  scaling  space  for  variables 

AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

nt  =  0 

do  12  i=l,ntot 

do  13  j  =  l ,  ntot 

if ( i .eq. j)then 

go  to  13 
else 

nt=nt+l 

xs(nt)=x(i)-x( j) 
endif 
13  continue 

12      continue 

nt  =  0 

do  14  i=l,ntot 

do  15  j  =  l , ntot 

if ( i .eq. j)then 

go  to  15 
else 


nt=nt+l 
ys(nt)=y(i )-y( j) 


endi  f 
15  continue 

14      continue 

DO  40  1=1,13 

COUNT( I )=0.0 
40      CONTINUE 

NTO=NTOT**2-NTOT 


AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

A 

*  This  section  of  the  program  calculates  the  ' z'    coordinate 

*  for  the  PDF.   The  'xmin'  and  'xmax'  values  are  the  beginning 

*  and  end  points  for  each  grid  in  the  x-y  plane. (  Similarly 

*  for  'ymin'  and  'ymax') 

A 
A 

DO  60  1=1,101 

DO  61  J=l,101 
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60      CONTINUE 

DO  50  1=1, NTO 

XMIN=-5.0 

XMAX=-4.9 

DO  51  K-1,101 

IF(XS( I ) .GE.XMIN.AND.XS( I ) . LT . XMAX ) THEN 
1XC  =  K 

GO  TO  5  2 
ELSE 

XHIN=XMIN+. 1 

XMAX=XMAX+. 1 
ENDIF 

51  CONTINUE 

52  YHAX=-5.0 

YMIN=-4 .9 

DO  53  K=l , 101 

IF( YS( I  ) . GE.XMIN.AND. YS( I ) . LT . YMAX ) THEN 

IYC=K 

GO  TO  54 
ELSE 

YMIN=YMIN+. 1 

YMAX-YMAX+.l 
ENDIF 

53  CONTINUE 

54  AMP( IXC, IYC)=AMP( IXC, IYC)+1 
50      CONTINUE 


DO  66  1=1, 101 

XX( I  )  =  I 

yy(i)-i 

66      CONTINUE 


CALL  SRFACE(AMP, 101,101,101 , XX , YY , WORK , S , 0 . 0 ) 
CALL  CLSGKS 


STOP 
END 
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